Model-informed health and socio-economic benefits of enhancing global equity and access to Covid-19 vaccines

We take a model-informed approach to the view that a global equitable access (GEA) to Covid-19 vaccines is the key to bring this pandemic to an end. We show that the equitable redistribution (proportional to population size) of the currently available vaccines is not sufficient to stop the pandemic, whereas a 60% increase in vaccine access (the global share of vaccinated people) would have allowed the current distribution to stop the pandemic in about a year of vaccination, saving millions of people in poor countries. We then investigate the interplay between access to vaccines and their distribution among rich and poor countries, showing that the access increase to stop the pandemic gets minimized at + 32% by the equitable distribution (− 36% in rich countries and + 60% in poor ones). To estimate the socio-economic benefits of a vaccination campaign with enhanced global equity and access (eGEA), we compare calibrated simulations of the current scenario with a hypothetical, vaccination-intensive scenario that assumes high rollouts (shown however by many rich and poor countries during the 2021–2022 vaccination campaign) and an improved equity from the current 2.5:1 to a 2:1 rich/poor-ratio of the population fractions vaccinated per day. Assuming that the corresponding + 130% of vaccine production is made possible by an Intellectual Property waiver, we show that the money saved on vaccines globally by the selected eGEA scenario overcomes the 5-year profit of the rights holders in the current situation. This justifies compensation mechanisms in exchange for the necessary licensing agreements. The good news is that the benefits of this eGEA scenario are still relevant, were we ready to implement it now.


S2.1 Epidemiological parameters in HI and MLI countries
We investigate if the Covid-19 evolution is similar in high-, middle-, and low-income countries, or if different epidemiological parameters should be adopted among the different groups of nations.Since we extend the model proposed by Gatto et al. [21] for the initial spread of Covid-19 in Italy (a rich country), we first test the hypothesis that the epidemiological parameters for middle and low-income countries are similar to those used for rich nations.To this end, we compare the strength and temporality of the Covid-19 outbreaks in the three groups of high-, middle-, and low-income countries.For each country (place of the outbreak), an outbreak is identified by a peak of the number of weekly new cases per million citizens.More precisely, the strength of the outbreak is a daily local maximum of the new confirmed cases of the last seven days (only considering maxima above the epidemiological threshold of 150 cases per million per week [57,58]) that is

Variable Definition
Si, S i susceptible (non infected, unprotected/partly protected) Ei, E i exposed (infected during latency, non infectious) Pi, P i presymptomatic (infectious, prior to possible symptoms onset) Ii unconfirmed symptomatic (infectious) I Q i quarantined symptomatic (infectious, isolated) Ai, A i untraced asymptomatic (infectious, not developing symptoms) A Q i quarantined asymptomatic (infectious, not developing symptoms, isolated) Hi hospitalized (infectious, isolated) Ri, R i unaware/aware resistant (non infected, fully protected) Di dead Table S1: State variables of the model.All variables are expressed as fractions of the world population.
The subscript i denotes the country of residence (i = 1 for HI countries; i = 2 for MLI ones).The primesuperscript denotes the classes of individuals who recently received a dose of vaccine and are developing resistance (except for the class R , that includes individuals aware of being resistant by vaccination or recovery from the disease).Classes I Q , A Q , and H include all subjects that are currently confirmed to be infected; they are all assumed to be isolated.from the previous peak (from the first case in the country for the first peak).for which the hypothesis of a full vaccine protection is questionable).A visual inspection highlights that data include outbreaks of moderate intensity for most low-income countries (red dots), whereas outbreaks are similar in middle-and high-income countries, as further highlighted by the box-plots reported in Fig.

S1(b)
. This is confirmed by statistical tests.A one-way ANOVA [59] (excluding outliers), confirms that there is a statistical difference between the outbreak average strength in the three groups (p-value of the no-difference hypothesis: < 10 −12 , p-value of the Kolmogorov-Smirnov test on residual normality: 0.09, pvalue of the Box test for residual homoscedasticity: 0.06).The multiple comparison marginal analysis [60] confirms that low-income countries have a significantly lower average with respect to the other two groups (95%-CI for the average of the three groups: low-income, [−221, 512] note the even negative lower bound; middle-income, [1797, 2996]; high-income, [2906,3596]).In the absence of relevant reasons to justify such a difference, and also excluding climatic factors (the same statistical difference is observed among countries at similar latitudes), we conclude that Covid-19 infection certification is not reliable (largely underestimated) in low-income countries.To further test the data similarity between high-and middle-income countries, we run the multi-variate MANOVA test [61] jointly on the strength and temporality of the outbreaks (p-value of the no-difference hypothesis: 0.12; p-values of the Kolmogorov-Smirnov test on data normality in each group: 0.65 and 0.24; p-values of the Box tests for homoscedasticity between the aggregated data (both groups) and the single groups: 0.56 for high-income countries; 0.33 for middle-income countries).We therefore conclude that there is no evidence for differences in the spread and evolution of Covid-19 in high-and middle-income countries, so that, discarding the unreliable data from low-income ones, we use the same epidemiological parameterization in the two macro groups.
Finally, because data on outbreaks do not explicitly account for mortality, we also compare the ratio of Covid-19 total deaths over the total number of confirmed cases across high-and middle-income countries (mortality data from low-income countries are not considered because unreliable, as well as data on contagions).The box-plots of the data, reported in Figure S2, suggest that no significant difference is present between the two groups.This is confirmed by an ANOVA test (p-value of the no-difference hypothesis: 0.17; p-values of the Kolmogorov-Smirnov test on data normality in each group: 0.32 and 0.24; p-values of the Box tests for homoscedasticity between the aggregated data (both groups) and the single groups: 0.14 for high-income countries; 0.10 for middle-income countries).These analyses allow us to conclude that neither mortality data show statistical evidences of difference in the two groups.
S2.2 Basic reproduction number R 0 , progress rates δ E , δ P , recovery rates γ I , γ A , γ H , and fraction of infected subjects developing symptoms σ We take from Gatto and coauthors [21] their estimates of the basic reproduction number R 0 = 3.6, of the progress rates δ E and δ P , of the recovery rates γ I , γ A , and γ H , and of the fraction σ of infected subjects developing symptoms (see Table S2).In particular, the latter parameter is rather uncertain, due to the untraceable nature of asymptomatic infections.Recently, a meta-analysis of several datasets distinguishing symptoms level among confirmed cases revealed a 40% fraction of totally asymptomatic subjects [62].The estimated fraction however misses the untraced subjects, so the asymptomatic fraction over all infections is expected to be larger (the value used by Gatto and coauthors is 1−σ = 0.75).Given the intrinsic uncertainty, we include parameter σ in the list of those for which a sensitivity analysis of our results is performed.
From R 0 , we derive the transmission rates β P 0 , β A 0 , and β I 0 in the absence of stringency measures.Specifically, we assume the same ratios between the transmission rates used by Gatto and coauthors and we use the same model formula for computing β P 0 from a given R 0 (see [21] Supporting Information), as it is computed at the disease-free equilibrium, where the models are equivalent.The resulting transmission rates are however slightly different because of different values of two other parameters: the isolation rate η (at home or in hospital) of symptomatic subjects and their mortality prior to isolation.The first has been increased, due to changes in people awareness and behavior from the early stage of the pandemic (the average time 1/η from symptom onset to isolation has been reduced from 4 to 1 day), while the second has been set to zero, neglecting those deaths (infected subjects can only die from class H).

S2.3 Fraction of hospitalized symptomatics 1 − ζ
To estimate the fraction of home treated (quarantined) symptomatic subjects ζ, we actually look at the complementary fraction 1 − ζ of hospitalized symptomatics.We use the weekly hospital admissions of Covid-19 patients and the weekly new Covid-19 cases, the latter corresponding in the model to the weekly inflow to the classes A Q , I Q , and H.However, parameter ζ only discriminates symptomatics in the model, so that, to remove asymptomatic cases, we use the estimate of 60% incidence of asymptomatic infections among the tested population reported in [62].Denoting the weekly inflow to a class with the 'in'-subscript, we can therefore link data and model by writing I Q in + H in = 0.4 (new weekly cases) and hence (new weekly hospitalized cases) 0.4(new weekly cases) .
In conclusion, we estimate 0.4(1 − ζ) as the least square linear regressor (with zero intercept) between the new weekly hospitalized cases and the new weekly cases.We used data from 21 high-income nations, for which both data series are available.The resulting regression is shown in Figure S3, yielding to the estimate ζ = 0.94.The determination coefficient R 2 = 0.23 is not so high, essentially because the fraction of hospitalized cases varied significantly between the first and successive infection waves and dropped substantially with the recent less severe virus variants (see dots below the regression line in the figure).Our model, however makes use of constant parameters, so that we keep the average fraction along the pandemic.

S2.4 Mortality rate α
In our model, infected subjects die only from the hospitalized class H (see eqs. (4)l-m).We therefore estimate the mortality rate α by looking at the ratio between daily deaths and hospital occupancy in the same day.
Data are available daily for 34 high-income nations.We apply a 7-day moving average and aggregate data over the considered countries to avoid weekly oscillations and attenuate other noise sources.The least square linear regression (with zero intercept) between these two series gives a slope of 0.019 (R 2 = 0.82).
However, data show a super-linear behavior, representing the fact that the mortality rate is higher when hospital capacity gets saturated (Fig. S4).We have therefore fitted a piecewise linear model, with a slope α

S2.5 Contact tracing rate a 0
As for the contact tracing rate a 0 at zero stringency measures, we qualitatively identify the ratio β P 1 /a 1 for the HI group on a small (but representative) dataset at the first level of stringency.This quantity is liked to a 0 by our modeling of the containment measures (see Sect.2).It indeed results β P 1 /a 1 = β P 0 /a 0 , independently of the imposed restrictions in group 1 (because the SL 1 -dependent factors in β P 1 and a 1 cancel out), from which we get a 0 = β P 0 /(β P 1 /a 1 ).Recall that β P 1 /a 1 is the ratio between two probabilities (because the daily contact rate included in both β P 0 and a 0 cancels out in the ratio): the probability of getting infected given a contact with a presymptomatic subject, and the one to be traced given a contact with a subject who is confirmed infected by the end of the day.It can be estimated as the average ratio between the number of secondary cases generated by a confirmed infection and the number of people that have been tested as a consequence of the primary case (e.g., if all contacts are traced, then the probability at the denominator is 1 and the ratio just gives the probability to get infected given the contact; otherwise it gives a larger value).
As a representative case, we found a local dataset on schools surveillance made public by Regione Veneto in Italy [56] (to the best of our knowledge this is the only data sample that allows the distinction of primary from secondary cases).Data report 4832 primary cases that caused 83660 tested people (students and school personnel) and 2793 secondary cases, in a stringency scenario imposing the use of masks and distancing in closed spaces.The resulting ratio is almost 30.The considered scenario is however one of the most severely traced in rich societies, so that we qualitatively set β P 1 /a 1 = 1/10, resulting in a 0 = 38.5.Because of the extreme locality of the dataset and given the uncertain role attributed to contact tracing on the long-term evolution of the epidemic [63,64,65,66], we include parameter a 0 in the list of those for which a sensitivity analysis is performed.

S2.6 Vaccination rates ρ 1 and ρ 2
Looking at Covid-19 vaccine administration, there is a significant difference between HI and MLI countries.In the status quo, i.e., the current scenario, vaccine production is the limiting factor, allowing the administration of a vaccine dose to about 0.5% and 0.2% of the population per day in HI and LMI countries, respectively, so that we set ρ 1 = 0.005 and ρ 2 = 0.002.Indeed, the available vaccination data show a clear difference between HI and MLI nations.In HI countries, the vaccination campaign proceeded with an average rate (over 2021) of 0.46 % of the population injected a dose per day (standard dev.0.18).We therefore set the rollout parameter ρ 1 to 0.005 in the status quo.
For MLI nations, we exclude China, because China decided to produce its own vaccine, calling itself out from the competition to acquire the vaccines on the private market.In the current scenario, China looks like an outlier, with very high vaccination rates in mid-2021 (see Fig. S5).Considering vaccination data from China would therefore raise the average vaccination rate in MLI nations, unrealistically showing the current inequity in access more globally equitable than it has been so far.Aggregating vaccination data in MLI nations (without China), see green line in Fig. S5, results in an average vaccination rate (over 2021) of 0.24 % of the aggregated population (standard dev.0.018).We therefore set the rollout parameter ρ 2 to 0.002 in the status quo.Note that the rollout parameters ρ 1 and ρ 2 are the control variables.Indeed, it represents the adopted vaccination campaign, which is decided and controlled by central governments.The rollout parameter defines both the access to vaccines ((N 1 ρ 1 + N 2 ρ 2 )/W = 0.16ρ 1 + 0.84ρ 2 ) and the inequity in access log 10 (ρ 1 /ρ 2 ), expressed in favor of HI countries.Clearly, these parameters highly influence the course of the Covid-19 pandemic.

High-income low-income low-income (without CHN) CHN
In the hypothetical scenarios, the control variable could assume whichever value.However, for the proposed eGEA scenario that we imagine feasible, we assume to reach 1% of the population injected per day (ρ 1 = 0.01, a peak value already realized with the current facilities), in HI countries.We assume that MLI countries can keep half of the rollout capacity of rich countries (ρ 2 = 0.005 a peak value already realized with the current facilities)).Note that this global daily access to vaccines ((N 1 ρ 1 + N 2 ρ 2 )/W = 0.0058) has been reached end even surpassed during the real vaccination campaign in specific days: the challenge is to constantly keep it..

S2.7 Acquisition (ν) and loss (ω) of virus resistance
We assume that vaccinated susceptibles develop full protection in about two weeks from each dose inoculation (1/ν = 14 days [36]).Regarding the duration of the protection, while there is evidence that antibody levels decline substantially within a few months from infection-asymptofaimatic or mildly symptomatic subjects with fewer antibodies [67,68]-vaccines are so far guaranteed to give protection for six months [37,38].At the qualitative level of this study, we optimistically assume that virus resistance wanes on average in 6 months after each dose or after recovery (ω = 1/180 day −1 , eqs.(4a,n,o)).

S2.8 Mobility probability C
Mobility between rich and poor countries is assumed low and symmetric (for simplicity, to avoid the consideration of migration flows).For all mobile classes Y ∈ M, we estimate the mobility probability as the fraction of EU citizens travelling (not for business) out of Europe (and not going to the USA) at least once in a year, resulting in approximately C = 0.0002 (source: ec.europa.eu/eurostat/cache/digpub/eumove).This is an underestimate of the probability C 12 for a citizen of group 1 (HI) to be in group 2 (MLI), as business trips are not counted and an average stay abroad of one day is considered.However, as we set a symmetric mobility (C Y 12 = C Y 21 = C), we purposely keep this underestimate and we keep it fixed, independently of the stringency index of the source and destination groups.

S2.9 ARX filter to reconstruct the SI
In this section, we report the details on the calibration of the autoregressive-exogenous (ARX) filter to reconstruct the stringency index (SI) from hospitalization data.Focusing on a collection of 19 high-income nations, we take daily hospitalization data of each nation (fraction of the country population in hospitals with confirmed Covid-19 infection, up to September 2021) as the filter input and the corresponding daily SI in the country as output (see the caption of Table S3 for the list of countries).We treated each national input-output series as an independent experiment and we used the Matlab Identification Toolbox to identify the ARX filter from the so-obtained multi-experiment dataset.
The model that minimizes the fit error reported by the Matlab Identification Toolbox makes use of 80 autoregressive daily steps and a weighted moving average on the last 10 days of hospitalization data, achieving a fit of 46% among the single national datasets.Indeed, two-three months is the typical duration of an infection wave, during which the SI is reasonably influenced by its own previous values, while 10 days of hospitalization data represents a typical monitoring window used by our policy-makers.The model is the following: where ŜI(t) is the SI reconstructed estimate at day t in the considered country (or average SI across countries), H(t) is the number of hospitalized subjects at the same day, and N is the population size.The coefficients α i and β i are reported in Table S3.Using the model on the aggregated data from the 19 considered nations (input signal: daily average of hospitalized population fractions; output signal: daily average SI), we obtain the result shown in the left panel of Fig. S6, where the standard deviation of the estimation error is 0.92 and the R 2 of the fitting is 0.99.
The auto-correlation and cross-correlation of the residuals are reported in the right panel, showing that the information carried by the input is fully exploited, as well as most of the information carried by the past values of the SI reconstruction (some information is left in the last 10 days of reconstruction).This is due to the intrinsic nature of the SI, since containment measures are not rediscussed daily, but rather on a weekly or longer basis.
We finally check if the SI evolves similarly in middle-income nations.We tested the ARX model (S1) on 20 middle-income nations, obtaining obtain the result shown in Fig. S7, a fit of 51%.The standard deviation of the estimation error is 4.24 and the R 2 of the fitting is 0.92.This result suggests a similar decision-making process in the two macro-groups of nations (assuming, in the absence of reliable data, that  S3 for the list of the 19 considered nations) and the reconstruction provided by filter (S1) fed by the average, over the same countries, of the population fractions of hospitalized subjects (blue dotted line).Right panel: Auto-correlation of the model residuals and their cross-correlation with the input signal.
low-income nations behave similarly to middle-income ones).It also provides a validation of the ARX filter itself.
We therefore implement two identical copies of filter (S1) in our model, one reconstructing the SI for HI countries and one for MLI ones.Taking the reconstructed value as the real one, we can formally remove the 'cap' from the filter output and add to our model the following two equations: We note, however, that hospitalization flows at the beginning of the simulation do not give rise to the highest stringency observed during the early stage of the pandemic.This is due to two concomitant factors: first, the hospitalized fraction of isolated symptomatics is estimated from data averaged across 2020-22 and is much lower than in the early stage (see SI Sect.S2.3); second, the initial fear for an unknown disease contributed to the highest stringency-with a worldwide nearly synchronous lockdown.To compensate for this discrepancy between the observed and simulated stringency, we force the highest stringency level SL i = 5 from day 46 (t = 45) for a period of 60 days, before and after which we use the stringency index provided by the filter.to appreciate substantial differences between this eGEA scenario and the status quo reported in Figure S8.
Most remarkably, the pandemic terminates in slightly more than a year of vaccination, while the pandemic stabilizes in an endemic regime with moderate waves in poor countries in the status quo. in the GEA scenario, averaged quantities are computed by averaging a zero value after the pandemic end and up to the 6-yr-horizon of the simulation).

S4 Sensitivity analysis
In what follows, we present the sensitivity analysis in absolute terms, separately for HI and MLI countries.
We see that, though contact tracing and mobility do not play a crucial role (up to doubling and halving their intensity), the incidence of symptomatic cases is more critical.Specifically, a higher incidence induces less overall cases in both scenarios, because symptomatic cases are quickly isolated, but more hospitalization and deaths.It also helps in stopping the pandemic earlier in the eGEA scenario, but only up to a small increase.If the incidence is doubled, for example, the pandemic is longer than it is for a 50% increase, probably because the reduction of cases limits the spread of the protection against the virus.The same trend is consequently observed in the required vaccine doses.
Finally, we present the sensitivity analysis with respect to the average duration ω − 1 (waning time) of vaccination-acquired or infection-acquired immunity.The effect of this parameter is crucial, giving a more severe pandemic the shorter is the duration is the waning time.Note however that the main feature of our proposed eGEA vaccination strategy, i.e., the fact that it is able to stop the pandemic in about a year of vaccination, is robust with respect to ω ∈ ( 1 360 , 1 90 ), meaning with respect to losing immunity on average in a range between 3 months and 1 year.
Overall, the sensitivity analysis confirms the robustness of our model calibration.The specific effects are however not all intuitive, essentially because of the intrinsic nonlinearity of the relations.This is indeed the added value of a model-informed approach.

S4.1 Vaccine efficacy
To relax the hypothesis of 100% vaccine efficacy, we need to consider that only a fraction e of susceptible subjects develop resistance against the virus, where the new parameter e is indeed the vaccine efficacy.
For modeling purposes, we assume that vaccinations on recovered resistant subjects (compartment R) have 100% efficacy, thus renewing their protection (flow R → R ), and that people on which the vaccine has been ineffective, face the same infection force of susceptibles.These latter subjects are unaware of their missing protection, so that they remain ineligible for a new vaccination for the average time ω − 1 of their presumed protection.Under such assumptions, we add a new compartment S for these latter subjects and consistently modify model ODEs (4) as follows: In the following table, the parameter e is reduced to analyze the sensitivity of our simulations of the status quo and eGEA scenarios with respect to the vaccine efficacy (e = 1 in the original model).Interestingly, note that the change in the pandemic length (for the eGEA scenario) has similar behavior to the one obtained in the sensitivity analysis of the protection waning rate ω for small parameter perturbations (compare the /1.03, /1.05, and /1.1 columns of this table with the ×1.03, ×1.05, and ×1.1 columns in the table of the sensitivity analysis for the protection waning rate ω).In other words, if less effective vaccines are used either because the protection time ω −1 is reduced (but all the vaccinated people become protected) or because their efficacy e is reduced (but protected people remain resistant for six months on average), the model simulations are analogous, both in terms of pandemic length and in the socio-economic indicators (i.e., case, deaths, vaccines doses, and average SI).Note however that, for intermediate perturbations (1.25 and 1.5), the model simulations foresee a big difference in the eGEA scenario: the pandemic stops perturbing the protection waning rate ω, but not perturbing the vaccine efficacy e.In other words, the sensitivity with respect to e is more critical with respect to the one with respect to ω.
This result can be intuitively explained.Indeed, by reducing the vaccine efficacy, we introduce in the model subjects who are susceptible to the infection while being not eligible for a new vaccination for an average time ω −1 .During this time, they fuel the epidemic, independently of the vaccination rate.Differently, by increasing the waning rate ω, people lose the acquired protection in a shorter time, but as soon as they lose it they are eligible for a new vaccination.The effect on the epidemic dynamics is therefore mitigated by an intense vaccination, as that proposed in our eGEA scenario.We now consider the effect of a reduced vaccine efficacy on all possible vaccination strategies, defined by our two control parameters: the global access to vaccines and inequity in access between HI and MLI countries.Figure S12 reports the total cases worldwide for the same reduced values of vaccine efficacy e considered so far.As the vaccine efficacy gets smaller, it is no longer true that an intense vaccination is able to stop the pandemic.This result holds true only for e > 65% (the first five panels of the figure, even if not visible in the fourth and the fifth panels).Under this condition, the minimum access granting the pandemic end is always obtained in a GEA scenario (for e = 80% the minimum access is 0.65%, while for e = 66% it is 2.1%, unfeasibly larger than the extreme 0.6% used in the figure).But the optimality of the GEA strategy at global scale remains confirmed independently of the vaccine efficacy.The GEA strategy minimizes the total cases worldwide (and, consequently, all the pandemic metrics at global scale) if a sufficient access to vaccines is granted (see the blue dotted line, that stabilizes around zero-inequity in the right part of each panel).

Figure S1 :
Figure S1: Covid-19 outbreaks in high-, middle-, and low-income countries.Panel (a) reports the average outbreak strength (measured in cases/milion/week, in log scale) versus the average outbreak temporality (measured as days from the previous peaks); green, blue, and red dots represent low-, middle-, and highincome outbreaks, respectively.Panel (b) reports the box plot of the average outbreak strength for the three groups.

Figure S2 :
FigureS2: Boxplots of the ratios between total deaths and cases among high-income (group 1 on the left) and middle-income (group 2) nations.

Figure S5 :
Figure S5: Time series of the vaccination rate (fraction of the population vaccinated per day) in the two groups of HI and MLI countries, along 2021.Data from China are separately plotted to show how China results as an outliers.

Figure S6 :
Figure S6: Left panel: Comparison between the averaged SI index available daily for high-income countries (black line, see TableS3for the list of the 19 considered nations) and the reconstruction provided by filter

Figure S10 :
Figure S10: Model simulation of the selected eGEA scenario (left/right: HI/MLI countries).The pandemic stops in about a year of vaccination.Note that the time horizon of five years of vaccination is long enough

ρ 1 /
Figure S12: Analysis of the total cases with respect to different vaccination strategies for different vaccine efficacy e.The dotted blue line shows the inequity that minimizes the cases at the specified value of access to vaccines.Vaccination strategies able to stop the pandemic before 6 years are encapsulated by the green line.

Table S2 :
Summary of the model parameters and their calibration.not exceeded in the previous and posterior 45 days.The temporality of the outbreak is the time (in days) ij 0.0002 Mobility probability Probability for an individual resident in group i in class X to be present in group j = i

Table S4 :
Pandemic metrics in the status quo and in the GEA scenario with access to vaccines as in the vaccination-intensive eGEA scenario considered in Tab. 1 (scaled to a world population W of 8 billion people;